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ABSTRACT 



Oh We describe computational tools that have been developed to simulate dynamical mass trans- 

fer in semi-detached, polytropic binaries that are initially executing synchronous rotation upon 
circular orbits. Initial equilibrium models are generated with a self-consistent field algorithm; 
models are then evolved in time with a parallel, explicit, Eulerian hydrodynamics code with no 
assumptions made about the symmetry of the system. Poisson's equation is solved along with 
the equations of ideal fluid mechanics to allow us to treat the nonlinear tidal distortion of the 
components in a fully self-consistent manner. We present results from several standard numerical 
experiments that have been conducted to assess the general viability and validity of our tools, 
and from benchmark simulations that follow the evolution of two detached systems through five 
full orbits (up to approximately 90 stellar dynamical times). These benchmark runs allow us to 
gauge the level of quantitative accuracy with which simulations of semi-detached systems can 
be performed using presently available computing resources. We find that we should be able to 
resolve mass transfer at levels M/M > few x 10~ 5 per orbit through approximately 20 orbits 
with each orbit taking about 30 hours of computing time on parallel computing platforms. 

Subject headings: binaries: close, hydrodynamics, methods: numerical 

1. Introduction which periodic or aperiodic variations in luminos- 
ity and spectral features can be explained by on- 
Over half of all stars m the sky are actually going mass-transfer events and instabilities in the 
multiple star systems and, of the binaries, about accretion flow. For example, long term stable mass 
half again are close enough to one another for transfer in which the accretor is either a white 
mass to be exchanged between the components dwarf; a neu tron star, or a black hole is widely rec- 
at some point in their evolution (Trimble 1983). ognized as the mechanism powering cataclysmic 
There is a subset of these close binary systems in variables (Warner 1995) and X-ray sources (Lewin, 



van Paradijs & van den Heuvel 1995). In Algol- 
type systems an evolved star transfers mass via 
Roche lobe overflow to a near main-sequence ac- 
cretor (Batten 1989; Vesper, Honeycutt & Hunt 
2001). Each of these systems evolves on a (secu- 
lar) time scale that is long compared to the orbital 
period of the system, with the mass-transfer rate 
determined by angular momentum losses from the 
binary, and thermal relaxation and nuclear evolu- 
tion of the donor star. In these systems, the frac- 
tion of the donor's mass that is transferred during 
one orbit is typically ~ 10~ 12 to 10~ 9 , many orders 
of magnitude less than what current numerical 3- 
D hydrocodes can resolve. All of the above sys- 
tems must have descended from binaries in which 
the accretor of today was initially the more mas- 
sive component who evolved first off the main- 
sequence. The mass transfer in these progenitor 
systems was in many instances dynamically un- 
stable, yielding mass-transfer rates many orders 
of magnitude above the currently observed rates 
and leading in some cases to a common envelope 
phase (see e.g. Warner 1995; Verbunt & van den 
Heuvel 1995; Nelson & Eggleton 2001) 

In addition, there is a wide class of binary star 
systems not presently undergoing mass-transfer 
for which the astrophysical scenarios that have 
been proposed to explain their origin or ultimate 
fate involve dynamical or thermal mass-transfer 
in a close binary, sometimes leading to a common 
envelope phase of evolution. Examples of such sys- 
tems are millisecond pulsars (Bhattacharya 1995), 
some central stars of planetary nebulae (Iben & 
Livio 1993), double degenerate white dwarf bi- 
naries, perhaps leading to supernovae of type la 
through a merger (Iben & Tutukov 1984), or 
subdwarf sdO, sdB stars (Iben 1990), and dou- 
ble neutron star binaries, perhaps yielding 7-ray 
bursts in a fireball when the neutron stars coa- 
lesce (Paczyhski 1986; Ruffcrt et al 1997; Mcszaros 
2001). The evolutionary scenarios that are drawn 
upon to explain the existence of some of these sys- 
tems call for events in which 10% or more of the 
donor's mass can be transferred during a single 
orbit. 

If we are to fully understand these rich classes 
of astrophysically interesting systems — their ori- 
gin, present evolutionary state, and ultimate fate 
— it seems clear that we will have to develop nu- 
merical algorithms that can accurately simulate 



mass-transfer events in binary systems under a 
wide range of physical conditions (for example, 
systems having a wide range of total masses, mass 
ratios, and ages) over both short and long evo- 
lutionary time scales. The astrophysics commu- 
nity as a whole is far from achieving this ulti- 
mate goal, but progress is being made as various 
groups are methodically tackling small pieces of 
this very large and imposing problem. Examples 
of recent progress in the numerical simulation of 
interacting binaries include two-dimensional simu- 
lation of mass transfer in Algol (Blondin, Richards 
& Malinowski 1995), three-dimensional evolutions 
of Roche Lobe overflow in LMC X-4 (Boroson 
et al 2001) and the accretion stream in /3 Lyrae 
(Bisikalo et al 2000), simulations of the common 
envelope phase and merger of a point-mass white 
dwarf with a red giant star (Sandquist et al 1998), 
neutron star binary and black hole-neutron star bi- 
nary mergers in the context of 7-ray bursts ( Janka 
et al 1999) and the dispersal of the secondary star's 
material in type la supernovae (Marietta, Burrows 
& Fryxell 2000). 

Building on our experience simulating the non- 
linear development of dynamical instabilities in 
self-gravitating systems — such as protostellar gas 
clouds (Cazes & Tohline 2000), stellar cores (New, 
Centrella & Tohline 2000), and young neutron 
stars (Lindblom, Tohline & Vallisneri 2001) - 
and on our experience in studying mass-transfer 
in close binaries through analytical and semi- 
analytical techniques (King et al 1997; Frank, 
King & Rainc 2001) we are developing a hydro- 
dynamical tool to study relatively rapid phases of 
mass-transfer in binary systems. Our immediate 
aim is to be able to follow the dynamical redistri- 
bution of material through ^10 — 30 orbits after 
the onset of a mass-transfer instability, in binary 
systems having a wide range of initial mass ra- 
tios with either star initially selected to be in con- 
tact with its Roche lobe and thereby become the 
"donor" . Our simulation tool treats both stars as 
self-gravitating fluids; they are embedded in the 
computational grid in such a way that their inter- 
nal structures are both fully resolved; and the sys- 
tem as a whole is evolved forward in time through 
an explicit integration of the standard fluid equa- 
tions, coupled with the Poisson equation so that 
the Newtonian gravitational field changes along 
with the mass distribution in a fully self-consistent 



way. Initially we will examine structures that can 
be well-represented by relatively simple barotropic 
(and adiabatic) equations of state, but this con- 
straint can easily be lifted in the future. We will 
be restricted to studies of relatively rapid phases of 
mass-transfer because we are integrating the equa- 
tions of motion forward in time via an explicit in- 
tegration scheme. 

While this simulation tool will not permit us 
to model stable flows with low mass-transfer rates 
- such as the observed flows in CVs and X-ray 
binaries — it should be capable of a wide range 
of astrophysical applications including: A deter- 
mination of the conditions required to become un- 
stable toward dynamical mass-transfer in all kinds 
of close binaries with normal and degenerate com- 
ponents; the ability to follow dynamical phases of 
mass transfer through to completion, which may 
mean a return to stability at a new system mass 
ratio, the formation of a massive disk or a common 
envelope with or without rapid mass loss from the 
system, or a merger of the two objects into one; 
and an investigation of the steady-state structure 
of secularly stable binaries. Through such inves- 
tigations we will be able to place on much firmer 
footing a variety of theoretical scenarios (as al- 
luded to above) that have been proposed to ex- 
plain the evolution and fate of close binary sys- 
tems. 

With the commissioning of gravitational wave 
interferometers such as TAMA, LIGO, and VIRGO, 
there has been a growing interest in understanding 
the detailed behavior of, especially, neutron star 
inspirals and mergers. As has been reviewed by 
Swesty, Wang and Calder (2000; hereafter SWC), 
a number of different groups have developed hy- 
drodynamical codes to simulate the late stages 
of inspiral and merger of such compact objects. 
Indeed, as has been described by New & Tohlinc 
(1997) and reviewed by SWC, an earlier version of 
our own simulation tool has been used to study the 
dynamical merger of equal- mass systems in which 
the stellar components were modeled with poly- 
tropic, white dwarf, and neutron star equations 
of state. Generally speaking, however, the last 
phase of a neutron star inspiral can be modeled 
with a hydrodynamical code that is less sensitive 
to initial conditions and more tolerant of errors in 
the algorithm that integrates the fluid equations 
forward in time than a hydrodynamical code that 



is designed to study more generic mass-transfer 
events in close binary systems. This is because 
general relativistic effects will necessarily drive a 
binary neutron star system to smaller separation, 
guaranteeing that the system will merge; and, 
even in the absence of relativistic effects, it ap- 
pears as though a tidal instability that disrupts 
one or both stars will be encountered before ei- 
ther fills its Roche lobe and encounters a classic 
mass-transfer instability (Lai, Rasio & Shapiro 
1994). 

Building on the work of New & Tohlinc (1997), 
we now have a simulation tool that can hydro- 
dynamically follow the orbital evolution of binary 
stars with high precision. In developing this tool 
we have made a number of improvements to the 
hydrodynamics algorithm that was used in this 
earlier work only to study the tidal merger prob- 
lem. We also have implemented a sclf-consistent- 
field algorithm that can construct very accurate 
initial equilibrium models of unequal-mass bina- 
ries in circular orbits, have paid special atten- 
tion to the manner in which initial models are in- 
troduced into the hydrodynamics code, and have 
taken full advantage of continuing improvements 
in high-performance computers. With this tool in 
hand, we should be able to accurately model the 
evolution of a much broader class of close binary 
systems, specifically, systems in which the compo- 
nents initially have unequal mass and/or radii and 
in which a mass-transfer instability rather than a 
tidal instability sets in. With the inclusion of ap- 
propriate relativistic corrections, this simulation 
tool should in principle also be able to simulate 
the merger of equal or unequal mass neutron star 
binaries, but our intention is not to focus so nar- 
rowly on this particular class of systems. 

In §2 of this paper, we collect results from theo- 
retical investigations of the linear stability of mass 
transfer in close binaries and discuss the approxi- 
mations that have been required to arrive at these 
results. In §3 we present the self-consistent field 
method we use for the construction of initial, equi- 
librium models. We then describe our implemen- 
tation of a parallel hydrodynamics code for the 
solution of the ideal fluid equations and Poisson's 
equation for an isolated mass distribution in §4. 
In §5 we compare results from the hydrodynamics 
code with known solutions for a set of test prob- 
lems and in §6 we present the results from the evo- 



lution of two benchmark detached binaries. These 
simulations demonstrate our ability to faithfully 
represent the forces acting on the fluid and allow 
us to estimate the mass transfer rate we will be 
able to resolve and the computational expense re- 
quired to evolve a given system through an inter- 
esting number of orbits. We conclude in §7 by 
summarizing the limits we have been able to at- 
tain at practical simulation resolutions and discuss 
the future application of the tool set to systems of 
interest. 

2. Theoretical Description of Dynamical 
Mass Transfer 

In this section we will argue that Roche lobe 
overflow in a binary system approximated by two 
polytropic components can result in mass trans- 
fer on a dynamical time scale for a certain range 
of polytropic indices. A spherical polytrope with 
uniform entropy in mechanical equilibrium obeys 
the following mass radius relation (c.f., Chan- 
drasekhar 1939), 



which implies, 
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Hence, the body will expand upon mass loss for 
polytropic indices satisfying 1 < n < 3. 

Consider Paczyhski's (1971) approximation for 
the effective Roche lobe radius R^ h of a donor 
star, taken to be the secondary, with mass M 2 in a 
point mass binary of total mass M and separation 
a, 
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From this, one obtains the following relation for 
the logarithmic change in the donor's Roche lobe 
radius, 

Rf h a 1M 2 1M 
a + 3Ma ~ 3M' 



i?? L 



(4) 



Upon eliminating the separation in favor of the 
system's orbital angular momentum J, one arrives 
at 

7?? L J 5M 2 2M Mi 

(5) 
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where M\ is the mass of the accreting star, taken 
to be the primary. If we further assume that the 
mass transfer is conservative with respect to the 
total mass and orbital angular momentum we de- 
duce that, 
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Comparing equation (2) with equation (6), the 
condition for stable mass transfer, R2 < i?^ L can 
be expressed as, 



& - £r > 0, 



(7) 



which, for a given polytropic index, implies a sta- 
ble mass ratio 
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(8) 



For a polytropic binary with n = 3/2 and mass ra- 
tio q > Stable = 2/3, mass transfer must occur on 
a dynamical time scale as the donor will readjust 
its structure within a few sound crossing times to 
its new mass. Note that if the donor is initially 
the less massive star (i.e. q < 1), the binary sep- 
aration is expected to steadily increase during the 
mass transfer event. But, if the donor is initially 
the more massive component (i.e. q > 1), conser- 
vation of orbital angular momentum implies that 
the separation must decrease and that the donor's 
Roche lobe radius will contract thus increasing the 
degree of overflow. The resulting mass transfer 
rate is expected in this case to be quite substan- 
tial. 

The dependence of the mass transfer rate on 
the degree of over-contact can be estimated from 
the product of the volume swept out by the flow 
near the inner Lagrange point, L\ 1 in unit time 
and the local value of the density. The cross sec- 
tion of the flow near L\ will scale as the square of 
the local sound speed, and the flow velocity is ap- 
proximately equal to the sound speed. The volume 
of material transferred in unit time then scales as 
the cube of the sound speed. The density near the 
edge of a spherical polytrope of index n, radius 
R 2 , mass M 2 , and polytropic constant k can be 
found by integrating the equation of hydrostatic 
equilibrium to obtain, 



p(r) 



GM 2 (R 2 - r) 



n(n + l) 
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(9) 



If we change variables to AR 2 = R 2 — r, the width 
of a spherical shell near the edge of the star, we 
obtain 

p(AR 2 )<x(AR 2 ) n . (10) 

The sound speed, c in turn obeys, 



so that 
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Taking the radius, r, to be the effective Roche lobe 
radius of the donor, AR 2 is the degree of over- 
contact. The mass transfer rate is then expected 
to scale with the degree of over-contact as 



dM 2 M 2 

T a 7 



AR 2 
~R~2~ 



(13) 



where P is the orbital period. This agrees with the 
calculation of Jedrzjec as presented in Paczyhski 
and Sicnkicwicz (1972). For a polytropic index of 
n = 3/2, eq. (13) indicates that the mass transfer 
rate will scale as the cube of the degree of over- 
contact. While the actual mass transfer rate ob- 
served in a fully self-consistent, three-dimensional 
evolution may differ substantially from the esti- 
mate given in (13), it nevertheless indicates that 
for unstable binaries mass-transfer events will 
evolve on a dynamical time scale once the donor 
reaches contact with its Roche lobe. 

All the results presented in this section have 
relied on a great many simplifying assumptions 
including the disregard of internal angular mo- 
mentum (spin) in each star, the use of the Roche 
model, neglecting the intrinsically nonspherical ge- 
ometry of the components, and assuming that the 
mass transfer event is, in fact, conservative. To 
proceed beyond this point one must deal with ex- 
tended distributions for the density and velocity in 
some approximation. We would argue further that 
it is advantageous to use a potential derived from 
the matter distribution in a self-consistent man- 
ner. With these additional complications, the task 
is well beyond the regime of analytical mechan- 
ics, but is tractable if we employ three-dimensional 
computational fluid dynamical techniques. To in- 
vestigate short time-scale mass transfer events nu- 
merically, we have developed a set of tools for both 
constructing equilibrium polytropic binaries and a 
hydrodynamics code to evolve systems of interest 
in time. These tools arc described below. 



3. Construction of Equilibrium Models 

The iterative method that we have used to 
generate equilibrium, polytropic binaries is very 
closely related to the self-consistent field (SCF) 
technique developed by Hachisu (1986; see also 
Hachisu, Eriguchi & Nomoto 1986). This tech- 
nique previously has been used to construct initial 
models of eguaZ-mass binary systems for dynami- 
cal studies of the binary merger problem (see New 
& Tohline 1997 and SWC). Here we employ a more 
generalized version of the technique to construct 
unequal-mass binaries. In the following discussion 
we use r to refer to an arbitrary point in space. 
The vector R is the cylindrical radius vector which 
can be expressed as R = ix + jiy. The axis of 
rotation is always taken to be parallel to, but not 
necessarily coincident with, the z-axis. 

The numerical results presented here are in a 
system of units where the gravitational constant, 
the radial extent of our numerical grid and the 
maximum density of one binary component are all 
taken to be unity. As we are treating polytropic 
models exclusively in the present work, these mod- 
els can be scaled to represent different physical 
systems by choosing a total system mass or or- 
bital separation for example. We would like to 
emphasize that the polytropic models could rep- 
resent binaries consisting of neutron stars, white 
dwarfs or normal stellar components. 

Assuming synchronous rotation so that the 
bodies are stationary in a corotating reference 
frame, the equations of hydrostatic equilibrium 
reduce to the following single vector equation, 



V I # + $- -ft 2 |R-R, 



= 0. 



(14) 



Here $ is the gravitational potential, £1 is the an- 
gular velocity of the reference frame in which the 
fluid is stationary, H is the enthalpy, and R CO m is 
the cylindrical radius vector of the system's center 
of mass so that |R — R CO m| is each fluid element's 
distance from the axis of rotation. For a poly- 
trope of index n, H is given to within an arbitrary 
constant by 



H 



— = (n + l) - = (n + l)npn , (15) 
P P 



where p and p are the pressure and density of the 
fluid, respectively. Equation (14) in turn implies 



that, 
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for some constant C, which in general will be dif- 
ferent for each binary component. Hereafter we 
denote these two integration constants as C\ and 
C 2 . 

Using equation (16) one can construct an iter- 
ative scheme as follows. An initial guess at the 
density field is constructed. Poisson's equation is 
solved to obtain the gravitational potential arising 
from the chosen mass distribution. This is, by far, 
the most computationally intensive part of the al- 
gorithm. For our work, we have chosen a cylindri- 
cal coordinate grid and utilized subroutines from 
the FISHPACK Fortran subroutine set for the 
solution of elliptic partial differential equations 
(Schwarztrauber & Sweet 1975; Schwarztrauber et 
al 2001) with the boundary potential being calcu- 
lated via a spherical harmonic expansion of the 
density field utilizing harmonic moments through 
£=10. 

With the gravitational potential in hand we can 
use algebraic relations at three boundary points 
where the density field is forced to vanish in order 
to set the integration constants C\ and C2, and 
the value of the angular velocity f2. The bound- 
ary points all lie along the line of centers. They 
correspond to the inner and outer boundary points 
for one star, and the inner boundary point for the 
companion as illustrated in Fig. 1. The values of 
the gravitational potential at the three boundary 
points, ta, tb and re are used to solve for the set 
of constants CI, C\ and C2 as follows: 
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C 2 = $(r c )-ir2 2 |R c -R com | 2 . 


(19) 



Equation (16) can then be used to construct the 
enthalpy throughout the computational domain 
and, from it, an improved density distribution can 
be constructed using the relation, 
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(20) 



where i = 1,2 labels the two stellar components. 
As Hachisu (1986) has explained, it is best to hold 
the values of p max ,i and p max ,2 fixed throughout 
the iteration cycles. 

The iteration cycle is then repeated using the 
improved density distribution until the relative 
change from iteration to iteration in C\, C2, CI, 
and -ffmax.i are all smaller than some prescribed 
convergence criterion, <5. For a grid resolution of 
128 radial points by 128 vertical points by 256 
points in azimuth, we typically use a tolerance 
5= 1 x 10~ 4 . 

Unfortunately the self-consistent field method 
does not allow one to specify physically meaning- 
ful parameters such as the binary mass ratio or 
separation a priori. Instead, as already described, 
it is best to specify the three boundary points and 
the maximum density for each body. Neverthe- 
less, the method described above remains, to our 
knowledge, the most effective means of generating 
fully self-consistent models of synchronously ro- 
tating, equilibrium binary systems with unequal 
masses and/or radii. 

We gauge the quality of a converged solution 
by the degree to which it satisfies the scalar virial 
equation. Specifically, we define the following di- 
mcnsionless virial error 
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(21) 



where the terms appearing in equation (21) are 
defined by the following integral quantities: 



K= X - J ,,vv<IV. 



W 
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(22) 



(23) 



(24) 



where v is the velocity field as measured in the 
inertial frame of reference. As applied in the SCF 
technique, the velocity is entirely due to the rota- 
tion of the frame, that is 



n= P dv, 



v = £1 x (R - R c 



(25) 



In Fig. 2 we plot density contours in the merid- 
ional plane for one contact binary system, three 
semi-detached systems, and two detached systems 



that we have constructed using the SCF technique. 
The more massive component always appears on 
the left of the plots. Figure 3 shows contours 
in the equatorial plane for the same six systems. 
The solid lines are at mass density levels of 10~ 5 , 
1CT 4 , 1CT 3 , 1(T 2 and 10 _1 , where the density 
has been normalized to the maximum density for 
each model, and the dashed line follows the self- 
consistently determined critical Roche surface for 
the system. The binaries all have a polytropic in- 
dex n = 3/2; other key parameters for these mod- 
els are listed in Table 1. Throughout this work 
we take the secondary component (denoted by a 
'2') to be the component closest to contact (clos- 
est to being the donor). The values of q = M 2 /Mi 
shown in Table 1 (as well as in Figs. 2 and 3) give 
the ratio of the mass of the secondary to the pri- 
mary; the stellar radii (Ri and R2) and Roche lobe 
radii (Rf L and i?^ L ) have all been normalized to 
the orbital separation. The stellar radii and Roche 
lobe radii are the radii of spheres that have a vol- 
ume equal to the star or critical Roche surface, 
respectively. For Model 4, the Roche lobe of the 
primary extends beyond the computational grid so 
the effective radius of its Roche lobe is only a lower 
limit. All six of these models were constructed on 
a cylindrical grid of 128 radial and vertical zones 
by 256 azimuthal zones. 

Table 2 lists the resulting virial error for the 
contact system (Model 1 from Table 1) con- 
structed on grids of differing resolutions. As the 
convergence criterion, 5, is decreased, the number 
of required iterations increases. For fixed resolu- 
tion, the overall quality of the solution does not 
significantly improve beyond some limiting value 
of 5, regardless of the number of iterations taken. 
As the resolution is increased, the virial error de- 
creases roughly in proportion to the square root 
of the number of grid points. 

Due to the symmetry of these initial models 
about the equatorial plane, we only calculate the 
models in the half space of z > 0. Assuming the 
line of centers coincides with the x axis, the tidal 
distortion of each star also is symmetric about 
the y — plane. Hence, further computational 
efficiency could be obtained with this technique 
by limiting the computational grid to only extend 
from to 7r in azimuth. To date, we have not 
enforced this additional symmetry constraint, al- 
though in practice the converged models display 




Fig. 1. — Position of the three boundary points in 
the equatorial plane for a SCF binary model. At 
the three boundary points (denoted by asterisks 
and labeled A, B, and C) the density is forced to 
vanish. The contours represent density levels for 
the converged model. 









Fig. 2. — Slice through the meridional plane for 
six example SCF binaries. Detailed parameters for 
these models are provided in Table 1. The solid 
contours are in the logarithm of the normalized 
density at levels of 1CT 5 , 1CT 4 , 1(T 3 , 1CT 2 , and 
10 -1 . The dashed curve traces the critical surface 
for the self-consistent Roche potential. The den- 
sity of the stellar components falls rapidly near the 
surface and on the scale used for these figures, the 
lowest two contours of the density coincide with 
one another. 



Table 1 
Initial Equilibrium Binary Models 



Model 


i 


max 
Pi 


Ri 


Rf L 


max 
Pi 


Ri 


R^ 


VE 


1 


1.0000 


1.00 


0.3720 


0.3723 


1.00 


0.3720 


0.3723 


1.5 x 10~ 4 


2 


1.2111 


1.00 


0.3056 


0.3580 


0.60 


0.3893 


0.3915 


1.4 x 10~ 4 


3 


0.4801 


1.20 


0.3727 


0.4401 


1.00 


0.3126 


0.3129 


3.4 x 10~ 4 


4 


0.1999 


1.00 


0.3817 


> 0.5194 


0.77 


0.2476 


0.2478 


2.8 x 10~ 4 


5 (EB) 


1.0000 


1.00 


0.2984 


0.3778 


1.00 


0.2984 


0.3778 


2.0 x 10~ 4 


6 (UB) 


0.8436 


1.20 


0.3180 


0.3919 


1.00 


0.3200 


0.3620 


2.2 x 10~ 4 



Table 2 

Convergence for SCF Method 

R z (j> 8 VE 

64 64 128 1.0 x 10~ 3 1.0 x 10~ 3 

1.0 x 10~ 4 6.0 x 10~ 4 

1.0 xlO" 5 5.5 x 10~ 4 

1.0 x 10~ 6 5.5 x 10~ 4 

128 128 256 1.0 x 10~ 3 6.9 x 10~ 4 

1.0 x 10~ 4 2.0 x 10~ 4 

1.0 x 10~ 5 1.5 x 10~ 4 

1.0 x 10~ 6 1.4 x 10~ 4 

256 256 512 1.0 x 10~ 3 6.3 x 10~ 4 

1.0 x 10~ 4 1.0 x 10~ 4 

1.0 x 10~ 5 5.2 x 10~ 5 

1.0 x 10~ 6 4.7 x 10~ 5 

1.0 x 10~ 7 4.7 x 10~ 5 



this symmetry. 

The SCF method is insensitive to the functional 
form for the initial guess of the density distribu- 
tion. For uniform spheres and spherically sym- 
metric Gaussian density distributions one can ob- 
tain the same converged model to machine accu- 
racy. We also note that we have found that more 
rapid convergence for models with soft equations 
of state (e.g., n > 3/2) can be achieved by using 
an even mixture of the current and previous po- 
tentials during the iteration. For more rigid equa- 
tions of state, where there is more mass at the 
boundary points and hence a greater coupling be- 
tween the solution near the boundary points and 
the global solution, such mixing is not necessary 
and the solution converges rapidly. 

4. Hydrodynamics Implementation 
4.1. Continuum Mechanics Formalism 

We have developed an explicit, conservative, 
finite-volume, Eulerian hydrodynamics code that 
is second-order accurate in both time and space 
to evolve the equilibrium binaries. The program 
is similar to the ZEUS code developed by Stone 
& Norman (1992). The integration scheme is de- 
signed to evolve five primary variables that are 
densities of conserved quantities: the mass den- 
sity, p, the angular momentum density, A, the ra- 
dial momentum density, S, the vertical momen- 
tum density, T, and an entropy tracer, r. The 
entropy tracer, 



t = (ep)-> , 



(26) 



where e is the internal energy per unit mass and 7 
is the selected ratio of specific heats of the gas. It 
is related to the entropy of the fluid through the 
relation, 

s = c p ln-, (27) 

P 

where c p is the specific heat at constant pressure. 
Using the entropy tracer in lieu of the internal 
energy per unit mass or the total energy density 
allows us to avoid the finite difference represen- 
tation of the divergence of the velocity field that 
must otherwise be used to express the work done 
by pressure on the fluid. 

For the evolutions presented in this paper we 
have set 7 = 1 + 1/n. We note, however, that by 



allowing the compressible fluid system to evolve 
with an adiabatic exponent that differs from this 
value, the stars will not be homentropic. For ex- 
ample, by selecting an appropriate value of 7 > 
1 + 1/n we can effectively model stars that are 
convectively stable and that obey a mass-radius 
relation quite different from the normal polytropic 
one specified by eq. (2), that is, different from 
£s = (1 — n)/(3 — n). By doing this, we expect 
to be able to closely approximate the mass-radius 
relationship of main sequence stars. A more de- 
tailed discussion of this idea is beyond the scope 
of this paper. 

The set of differential equations that we solve is 
based on the conservation laws for these five con- 
served densities. Mass conservation is governed by 
the continuity equation, 



f + V. ( pv)=0, 



(28) 



where v is the velocity field. The velocity vector 
is expressed in terms of its components in a cylin- 
drical coordinate system as v = !iejj+u e^+w^. 
The three components of Euler's equation govern 
changes in the momentum densities. We express 
these equations in a frame of reference rotating 
with a constant angular velocity, fi, as follows: 



f + V.(5v, 



p^nr + tt^ + ^ o> ( 29 ) 
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dA 
~dt 
where, 



d<J> cff 
V • (Aw) = -p^r— - 2nSR, (31) 
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(32) 



The second and third terms appearing on the 
right-hand side of eq. (29) represent the curvature 
of cylindrical coordinates and the radial compo- 
nent of the Coriolis force, respectively. Likewise, 
the last term appearing on the right-hand side of 
eq. (31) represents the azimuthal component of the 
Coriolis force. 

From the first law of thermodynamics we know 
that in the most general case, the entropy tracer 
obeys the expression, 



dr t ds 

dt c p dt 



(33) 



Here we will be considering only adiabatic flows, in 
which case ds/dt = 0, so the entropy tracer obeys 
an advection equation of precisely the same form 
as the continuity equation, namely, 



Or 

dt 



V • (rv) = 0. 



(34) 



Even though we are performing adiabatic evolu- 
tions we can not simply use an adiabatic equation 
of state (p = k/9 7 ) and disregard the first law of 
thermodynamics because the polytropic constant 
is, in general, different for each binary component. 
Finally, we solve Poisson's equation once every 
integration timestep in order to calculate the force 
of gravity arising from the instantaneous mass dis- 
tribution, 

V 2 $ = AirGp; (35) 

and we use the ideal gas law as the equation of 
state to close the system of equations, namely, 



P = (7 - 1) t1 = (7 - 1) Pe- 



(36) 



It may be argued that our treatment of the ther- 
modynamics of the system as the purely adiabatic 
flow of an ideal fluid is overly simplified. How- 
ever, we believe that the self-consistent treatment 
of both binary components in the presence of the 
full nonlinear tidal forces is sufficiently complex 
and novel to warrant the use of a simple equation 
of state at the present time. This will allow us 
to establish the qualitative behavior of systems in 
this limiting case before additional complications 
leading to nonadiabatic heat transport are intro- 
duced into the simulations. 

4.2. Finite Volume Representation 

Before proceeding with the discussion of the hy- 
drodynamics algorithm that we have implemented 
to solve the equations presented in §4.1 we first 
describe the discretization that has been used to 
represent the exact partial differential equations 
when they are expressed as approximate algebraic 
relations between discrete points in the computa- 
tional grid. As in the ZEUS code, all scalar vari- 
ables and the diagonal components of tensors are 
defined at cell centers. The components of vectors 
are defined at the corresponding faces of the cell. 
A volume element and the relative positions of the 
variables within each cell is illustrated in Fig. 4. 
The cell extends from Ri to Ri+i in radius, from 



Zj to Zj+i in the vertical coordinate, and from <f>k 
to <f>k+i in the azimuthal coordinate. We repre- 
sent the staggered variables in the computational 
mesh with a half-index notation; the coordinates 
of the center of a grid cell are given by i? i+1 / 2 , 
z j+i/2i 0fe+i/2i f° r example. A complete listing of 
the variables and their centering is given in Table 
3. 

4.3. Treatment of Advection Terms 

Through the method of operator splitting, one 
can construct a numerical scheme that groups 
terms of the same physical character together. 
Again, following along the lines of the ZEUS code 
we implement a splitting scheme that separates 
updates of the fluid state due to Eulerian trans- 
port (advection) from updates due to the source 
terms. In this section we describe our treatment 
of the advection terms. 

Given the density A of any conserved quantity 
A that satisfies a generic conservation law of the 
form, 

|_ + V-(Av)=0, (37) 

we can replace the differential equation (37) with 
an equivalent integral equation, 



0_ 



f \dV = - [ V-(\v)dV = - [ AvdS. 

Jv Jv Js(v) 

(38) 
Equation (38) must hold for any volume. In par- 
ticular, it must hold for every volume element 
within the computational grid. The exact integral 
relation is then expressible in the following finite 
volume form for each grid cell: 

\ (n+advection) \t n ) 



1 



At 



_^A*vASi, (39) 



i=\ 



where the summation is over all six faces on the 
surface of the three-dimensional cell. The surface 
elements, ASi, are naturally face-centered with 
respect to the control volume in question, so av- 
erages must be taken to obtain the advection ve- 
locity components necessary to perform the dot 
product for the momentum densities as shown in 
eq. (39). We use second-order accurate, linear av- 
erages to construct the advection velocities in this 
case. The amount of A advected through each face 
is given by an upwind biased, linear interpolation 
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Table 3 
Hydrodynamic Variables and their Centering 



Centering 


Variable 


Description 


A 


Ri+l/2 


Cylindrical Radius Coordinate 




Zj+1/2 


Vertical Coordinate 




4>k+l/2 


Azimuthal Coordinate 




Pi+l/2,j+l/2,k+l/2 


Mass Density 




T i+l/2,j + l/2,k+l/2 


Entropy Tracer 




Pi+l/2,j + l/2,k+l/2 


Pressure 




Hi+l/2,j + l/2,k+l/2 


Enthalpy 




^i+l/2.j + l/2,k+l/2 


Gravitational Potential 




Qll 

^i+l/2,j+l/2,k+l/2 


Diagonal Components of Artificial Viscosity 


B 


Si,j+l/2,k+l/2 


Radial Momentum Density 




u iJ + l/2,k+l/2 


Radial Velocity 


C 


Ti+l/2,j,k+l/2 


Vertical Momentum Density 




w i+l/2,j,k+l/2 


Vertical Velocity 


D 


A-i+l/2,j + l/2,k 


Angular Momentum Density 




v i+l/2,j + l/2,k 


Azimuthal Velocity 





Fig. 3. — Slice through the equatorial plane for 
the same six binaries shown in Fig. 2. The solid 
contours are in the logarithm of the normalized 
density at levels of 1(T 5 , 1(T 4 , 1(T 3 , 1CT 2 and 
10 _1 . The dashed curve traces the critical surface 
for the self-consistent Roche potential. 



of the distribution of A to give A* as described by 
van Leer (1979). By construction, the amount of 
A that is transported out of one cell immediately 
flows into the neighboring cell; thus ensuring the 
conservative nature of the advection scheme. 

Unlike the ZEUS code, we do not use operator 
splitting along the three separate dimensions dur- 
ing the advection step. Instead, we perform the 
updates due to advection in all three dimensions 
simultaneously. We thus avoid concerns about 
bias that may be introduced by using an unsym- 
metrized ordering of the advection sweeps. A dis- 
cussion of how we obtain second-order accuracy in 
time for the advection step through time centering 
of the terms appearing in eq. (39) is presented in 
§4.6. 

Our advection scheme automatically reverts to 
a first-order accurate (upwind) scheme at local ex- 
trema in the primary fluid variables. In addition, 
it is necessary to introduce an artificial viscosity 
to stabilize the scheme in the presence of shocks. 
The artificial viscosity prescription we have imple- 
mented is detailed in 84.5. 
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4.4. Treatment of Source Terms 




c j+l 



Fig. 4. — Volume clement for a cell-centered quan- 
tity (denned at the open circle labeled A) in our 
cylindrical coordinate mesh. Radial, vertical and 
azimuthal face-centered quantities are defined at 
points B, C, and D, respectively. Table 3 lists all 
of the variables used in the hydrodynamics code 
along with their centering according to this dia- 
gram. 



The Lagrangian source terms for the momenta 
that are shown on the right-hand sides of eqs. (29)- 
(31) arise from the forces of pressure and gravity, 
as well as from the differentiation of the curvilin- 
ear basis vectors and the rotation of the reference 
frame. We have found it advantageous to com- 
bine the pressure gradient with the gradient of the 
gravitational potential, which results in a gradi- 
ent of the sum of H and <f>. Since the centrifugal 
force can also be expressed as the gradient of a 
potential, it is included as well to form an effec- 
tive potential as defined in eq. (32). As explained 
in §3, our initial models have the property that 
$ eff = constant everywhere, hence to reasonably 
high precision V$ cff = throughout both stars 
initially. 

The expressions we have used for the source 
term updates of the momentum densities are given 
here by expressions (40)- (42). As with the ad- 
vection step, we do not use an operator splitting 
technique to evaluate the source terms along the 
three separate coordinate dimensions; instead, at 
each cell location, all updates due to Lagrangian 
source terms are performed simultaneously. 



-»(n+source) 



s 



(n) 



i,j+l/2,k+l/2 ij + l/2,/c+l/2 



At 



Pi,j+l/2,k+l/2 



AR 



(40) 



^i+l/2,j+l/2,k+l/2 ^i-l/2,j + l/2,k+l/2 



A 



i,j+l/2,k+l/2 
Pi,j+l/2,fc+l/2-Rj 



2flA 



i,j + l/2,k+l/2 _ 



T. 



(n+sourcc) 



T. 



(„) 



i+l/2,j,fe+l/2 «+l/2j,fc+l/2 



At 



(41) 



Pi+l/2,j,k+l/2 

~Az 



<f> cff — <b cS 

^i+l/2,j + l/2,k+l/2 ^i+l/2,j-l/2,k+l/2 



.(n+sourcc) A^ 

i i+l/2J + l/2,k ~ i+l/2,j+l/2,fe 



At 



(42) 



Pi+l/2,j + l/2,k 
A(j) 



<b cS — <b cS 

^i+l/2,j + l/2,k+l/2 ^i+l/2,j + l/2,k-l/2 



— 2lXSj + i/2,j+l/2,fc-Ri+l/2- 

Note that a caret identifies a variable whose value 
has been interpolated to a spatial location that 
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is different from the variable's primary definition 
point as shown in Fig. 4. These variables are given 
by volume-weighted averages as follows: 



A 



1 
i,j+i/2,fc+i/2 - ^r 

|_(^i+l/2,j+l/2,fe + ■^■i+l/2,j+l/2,k+l) 
1 



(43) 



Ri 



-AR 



+ {■^-i-l/2,j+l/2,k + ^i-l/2,j+l/2,fc+l) 



Ri - ^ARj 



S. 



1 



i+l/2,j + l/2.k 



4iL. 



-1/2 



(44) 



[(^»+l,j+l/2,fe+l/2 + Si+lJ+l/2k-l/2) 

1 



Ri- 



■1/2 



:AR,, 



-1/2 



+ \Si,j+l/2,k+l/2 + S'jj + l/2,fc-l/2) 
Ri+1/2 - J^R 



1 

2i?" 



Pi,j+l/2,k+l/2 



Pi+l/2,j+l/2,k+l/2 [ Ri + T^-R 



(45) 



+ Pi-l/2,j+l/2,k+l/2 I -R-. 



-AR 



Pi+l/2,j,k+l/2 — X (ft+l/2j + l/2,fe+l/2 (46) 
+ Pi+l/2,j-l/2,fe+l/2) j 

Pi+l/2,j+l/2,fc = 2 (P»+l/2,j+l/2,fe+l/2 (47) 
+ Pi+l/2,j+l/2,k-l/2) ■ 

4.5. Artificial Viscosity 

To stabilize the scheme in the presence of 
shocks, we employ a planar, von Neumann artifi- 
cial viscosity that is active only for zones that are 
undergoing compression. (See Stone & Norman 
1992 or Bowers & Wilson 1991, pg. 311 for more 
detailed discussions of artificial viscosity in Eule- 
rian hydrodynamics.) The momentum densities 
are updated from the following finite-difference 
equations, 



S 



(n+viscosity) 
,j+l/2,k+l/2 



»,j+l/2,fe+l/2 



At 



1 

AR 



(48) 



f)RR 
\yi+l/2,j+l/2,k+l/2 

rp(n+viscosity) rp(n) 
i+l/2,j,k+l/2 i+1/ 


n RR 
{ °i!i-l/2,j+l/2,k+l/2 

2,j,k+l/2 1 


At 

[Qi+l/2,j + l/2,k+l/2 ~ Qi+l/2,j- 

\ (n+viscosity) * (n) 
i+l/2j + l/2,fc ~~ i+l/2,j+l/2.fc 


Az 

-l/2,fe+l/2 
1 




At 




Acj) 



(49) 



(50) 



(q 



i+l/2,j + l/2,k+l/2 



-Q 



i+l/2,j+l/2,k-l/2 J ' 



where the diagonal components of the artificial 
viscosity are given by, 

n RR 

{ *i+l/2,j+l/2,k+l/2 



Qi+i 



+ l/2,j+l/2,k+l/2 



Q 



i+l/2,j + l/2,k+l/2 



VPi+1/2 


j + l/2,k+l/2 




(51) 


\ u i+l,j+l/2,k+l/2 — 


u i,j+l/2,k+l/2) j 


VPi+1/2 


j + l/2,k+l/2 




(52) 


(wi+1/2 


.j + l.k+1/2 - 


w i+l/2, 


j,k+l/2) j 


VPi+l/2 


j + l/2,k+l/2 




(53) 



( w i+l/2,j + l/2,fe+l — v i+l/2,j + l/2,k) 



if the velocity difference is negative; otherwise the 
components of Q are zero. Note that we neglect 
the shear components of viscosity. The factor v is 
a parameter that roughly dictates the number of 
zones across which shock structures will be spread. 
A value of v = 2 is typically sufficient. In keeping 
with our overall adiabatic treatment of the flow 
(see §4.1), we neglect the generation of entropy by 
shock compression. 

In a binary system that is undergoing mass 
transfer, the accretion stream will necessarily un- 
dergo a shock transition as it is decelerated upon 
impact with the accreting star, or when it inter- 
sects itself if the stream has sufficient angular mo- 
mentum to orbit the accretor. In addition, even 
for a detached binary simulation there will be weak 
standing shock fronts (as viewed in the corotating 
frame of reference) at or near the surface of the 
stars arising from the rapid deceleration of mate- 
rial falling onto the stars. We have found that even 
these weak shocks can have a noticeable impact 
on the quality of the solution in long time evolu- 
tions of detached systems unless artificial viscosity 
is used to damp the resulting oscillations. 

4.6. Time Centering 

The timestep cycle is split between the applica- 
tion of source, advection and viscosity operators. 
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First, the source terms are applied for one half of 
a timestep. Next, all updates due to advection 
are performed for a full timestep and the viscosity 
updates are applied to the momentum densities. 
Finally, the second half of the source operators 
are applied. The source and advection steps are 
thereby staggered in time when viewed over sev- 
eral iteration cycles for a constant value of the 
timestep. 

The advection is time-centered by first perform- 
ing half a timestep of fictitious advection in order 
to obtain "time-centered" velocities for construct- 
ing the face-centered advection velocity compo- 
nents that appear in eq. (39). The full timestep 
of advection is then performed. The components 
of the viscosity tensor are constructed from the 
velocity and density estimates at the midpoint of 
the timestep as well. 

Since the momentum densities themselves also 
appear in the source terms of eqs. (29) and (31), 
similar care must be taken with their centering in 
time. The source operators are applied in a ficti- 
tious source step to obtain the angular and radial 
momentum densities at a point half a timestep in 
the future. These values are then used to update 
the momentum densities through a full timestep. 
As the timestep value changes from iteration to 
iteration, this algorithm for time centering the 
source terms is not formally accurate to second 
order. However, in real computations the char- 
acter of the flow and, hence, the maximal signal 
velocity do not change rapidly over the course of a 
timestep cycle so that one may expect the result- 
ing inaccuracies in the time centering of the source 
terms to be small. The other terms that appear in 
the source operators, including the gravitational 
potential, are all calculated at the approximate 
midpoint in time between the source steps. 

4.7. Timestep Formulation and Boundary 
Conditions 

Since we explicitly integrate the fluid equations 
in time, the timestep is limited in size by the famil- 
iar Courant-Fricdrichs-Lewy (CFL) stability cri- 
terion which ensures the time increment is small 
enough so that no characteristic can cross a cell in 
a single timestep. Specifically, 



where c is the speed of sound. In practice we limit 
the timestep to a half the CFL time. Since we have 
introduced the diffusion terms associated with ar- 
tificial viscosity, the timestep must also satisfy the 
condition, (see p. 270 of Bowers & Wilson 1991), 



1 • 
At < —mm 

~ 4 



pAR pAZ P RA4> 



)RR' 



)ZZ ' 



1/2 



(55) 



At = 



AR AZ RA<j) 



(54) 



The boundary conditions for the fluid variables 
at the external boundaries are to allow the fluid 
to flow freely off the grid but to not allow material 
to flow back from the outermost layer of bound- 
ary cells. The central annulus of cells that has an 
inner radius at the coordinate axis is treated as a 
single azimuthally averaged cell for each layer in 
the vertical direction. 

4.8. Parallelization of Hydrodynamics Al- 
gorithm 

As it is our intention to perform high resolu- 
tion simulations, it is imperative that the work 
load within the simulation be distributed amongst 
many processors so that the simulations may be 
conducted in a reasonable amount of time and not 
exceed the available memory of a single node. The 
fluid dynamics equations, being hyperbolic partial 
differential equations, are ideally suited to a simple 
domain decomposition or single program multiple 
data (SPMD) parallelization model. Each com- 
putational task performs the same operations on 
their own block of the global data arrays with com- 
munication only being necessary between nearest 
neighbor tasks that share a boundary of ghost 
zones that is one-cell thick (this ghost zone thick- 
ness is dictated by the order of our advection and 
finite-difference operators). We have written the 
program in Fortran90 with explicit message pass- 
ing being performed with MPI (Message Passing 
Interface) subroutine calls. The resulting parallel 
code performance scales linearly with the number 
of processors for 4 to 128 processors on the Cray 
T3E. Similar behavior is also seen on the IBM SP 
platform. 

4.9. Solution of Poisson's Equation 

We are seeking to solve Poisson's equation for 
an isolated distribution of mass. The correct 
boundary condition in this instance is that the 
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potential goes to zero at infinity. As we only con- 
struct the solution on a finite domain we must 
specify the gravitational potential (or its gradi- 
ent) on some boundary that encloses all the mass 
in the simulation. We construct the boundary po- 
tential using a novel technique based on a com- 
pact representation of the cylindrical Greens func- 
tion in terms of half-integer degree Legendre func- 
tions of the second kind as described by Cohl & 
Tohline (1999). The boundary potential is then 
simply given by the convolution of the appropri- 
ate Greens function with the density distribution. 
This method is capable of generating the exact so- 
lution for a discretized mass distribution and has 
the attractive feature that it can be applied to 
very flattened bodies without suffering penalties 
in either performance or accuracy. 

In order to obtain the interior solution for the 
gravitational potential, Poisson's equation is first 
Fourier transformed in the azimuthal direction, 
then the resulting set of two-dimensional partial 
differential equations (Hclmholtz equations) for 
the decoupled Fourier amplitudes are solved us- 
ing an alternating direction implicit (ADI) scheme 
(c.f., Pcaceman & Rachford 1955; Black & Boden- 
hcimcr 1975). The solution is then transformed 
back to real space. 

The solution of Poisson's equation requires spe- 
cial care in the context of parallel computing be- 
cause the solution necessarily involves global com- 
munication as the character of the underlying 
physical law is action at a distance. The algo- 
rithms we have implemented for computing the 
gravitational potential are well suited to a cylin- 
drical geometry and very efficient in a distributed 
computing environment. Parallel communications 
are used to transpose the data so that all the 
data in a given dimension are in local memory at 
one time. When operations are to be performed 
along a different dimension, the data are trans- 
posed again. This allows us to send a relatively few 
number of large messages. Further details regard- 
ing our solution of Poisson's equation in a parallel 
computing environment can be found in Cohl, Sun 
& Tohline (1997). 

5. Test Cases 

Here we present results from three different 
types of tests that we have used to evaluate and 



quantify the accuracy of our computational tools. 
In all tests we compare a known, although not 
necessarily analytical, solution with the calculated 
numerical solution. 

5.1. Riemann Shock Tube Test 

As a check of the stability of our code in the 
presence of, ideally, discontinuous jumps in the 
fluid variables we have solved Sod's shock tube 
problem (Sod 1978) with the initial discontinuity 
lying along a plane of constant z — zq. Sod's 
shock tube problem presents a useful hydrody- 
namic test because the solution is known analyti- 
cally and contains the three simple waves that can 
occur in ideal fluid flow. Of these simple waves, it 
is the shock wave that concerns us most. Our goal 
is not to resolve the details of the shock structure 
but rather to ensure that our algorithm is well be- 
haved (numerically stable and yields an acceptable 
solution) in the presence of shocks. 

The initial conditions for Sod's shock tube 
problem are that the velocity, v, is zero every- 
where; for z < zo, the pressure, density and in- 
ternal energy per unit mass take on the values 
Pi = 1.0, pi = 1.0, ei = 2.5; for z > z , P u = 0.1, 
p u = 0.125, e u = 2.0. The fluid flow is character- 
ized by an adiabatic exponent, 7 = 1.4. 

The computed solution for the vertical veloc- 
ity, w, pressure, p, mass density p, and the quan- 
tity, t J p (which is proportional to the polytropic 
constant and, hence, the entropy; see eq. 27) is 
plotted along with the analytical solution at time 
t = 0.247 in Fig. 5. The computed points are 
not average values but are instead the values for 
a random column of cells at constant radius and 
azimuth within the three-dimensional grid. The 
calculation was performed with a coefficient for 
the artificial viscosity of v = 2.0 and with 130 ver- 
tical zones. The initial discontinuity was placed 
at zq = —0.1 and the grid extended from —1/2 to 
1/2 in the vertical direction. 

The results from this simulation compare favor- 
ably to the results produced by other second-order 
accurate, Eulerian hydrodynamics programs with 
artificial viscosity (c.f., Hawley, Wlison & Smarr 
1984; Stone & Norman 1992; Lufkin & Hawley 
1993). The shock front is spread out over ap- 
proximately three zones, and there is no indica- 
tion of numerical instability in the solution for 
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the shocked gas. The contact discontinuity is like- 
wise spread out over about three zones due to the 
numerical diffusion inherent in a second-order ac- 
curate Eulerian scheme. There is some disagree- 
ment between the computed and analytical so- 
lution at the tail of the rarefaction wave. This 
phenomenon has been investigated by Norman & 
Winkler (1983) and results from an inconsistent 
representation of the analytic viscous equations 
in finite difference form. Finally, immediately be- 
hind the shock, the analytical solution for t j p dis- 
agrees slightly with the computed solution for the 
shocked gas. This is because the analytical so- 
lution includes a small increase in the entropy of 
the fluid that passes through the shock whereas, as 
discussed in §4.5, we have elected to treat all of the 
flow as through it is adiabatic, that is, by setting 
ds/dt = in cq. (33). [Agreement with the ana- 
lytical result could readily have been achieved by 
constructing an appropriate expression for ds/dt 
in terms of the artificial viscosity in order to ac- 
count for dissipation in the shock, as shown for 
example in eqs. (11) and (37) of Stone & Norman 
1992.] In effect, we have assumed that each fluid 
element that passes through the shock is immedi- 
ately able to cool back down to a temperature that 
places it back on its original pre-shock adiabat. 

We note that we have used the gradient of the 
pressure as opposed to the density times the gradi- 
ent of the enthalpy for the solution of Sod's prob- 
lem. Due to the pathological nature of the discon- 
tinuous initial conditions, a correct solution can- 
not be obtained if the enthalpy term is used with 
our chosen centering of the fluid variables. 

5.2. Test of Poisson Solver 

Cohl & Tohline (1999) have published exhaus- 
tive tests showing the accuracy with which we are 
able to evaluate the gravitational potential on the 
boundary of our cylindrical coordinate grid. In 
order to ascertain the accuracy with which we 
are able to determine the force of gravity arising 
from the fluid everywhere inside the grid, we have 
calculated the potential and its derivatives for a 
uniform-density sphere of radius i?* and density 
p*, centered at an arbitrary position on the grid, 



tq. The analytical potential is 



$(r) 



-27rGp* [Rl - — d < R* 

4. k J (56) 

-— Gp*— ^- d>R*, 

3 d 



where d = |r — i"o I - The corresponding derivatives 
appearing in the gravitational force are: 

9$ _ AirGp* 
dR ~ ~ 



[(R cos (f> — xq) cos <fr 



+ (Rsmcf) — j/o) sin 

<9$ _ 4ttGp 

~dz~ ~ ~ 

<9$ A-KGp 



(z - zo) , (57) 

[-R (R COS (j) — Xq) Sin (j) 



R (R sin (f> — 2/0 ) cos <j)\ , 



for d < i?», and 

<9$ AirGp* Rl 



dR 

<9$ 
dz 



3 d? 
(Rsincj)- y )sm<f)} , 

AirGp* R 3 



[(R COS (j) — Xq) COS <fr 



3 d 3 
4ttGp* Rl 

~ ~ IF 



(z - z ) , (58) 

[-R (R cos cf> — xq) sin (j) 



+ R (R sin (f> — yo) cos 4>] , 

for d > i?». 

In Table 4 we present the average relative error 
in the potential and the average absolute error in 
the three derivatives for a uniform density sphere 
(p* = 1) of radius _R» = 1/3 placed at the origin 
and at xo = 1/2 for a representative set of grid res- 
olutions. The grid extends from to 1 in radius 
and from —1/2 to 1/2 in the vertical direction. 
Similarly, in Table 5 we present the maximum er- 
rors for the same quantities. 

The region near the surface of the sphere con- 
tains the largest errors in the potential solution 
(c.f., Stone & Norman 1992). At the surface, the 
density falls discontinuously to zero and the slope 
of the solution changes abruptly. When placed at 
the origin, this high error region is resolved by a 
larger number of smaller- volume cells than when 
the sphere is placed off-axis in the grid. This re- 
sults in a worse average error for the potential and 
its radial and vertical derivatives for the axisym- 
metric solution despite the fact that the maximal 
errors are generally smaller in this instance. 
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Table 4 
Average Error for Gravitational Potential and Force 



Origin 


R 


z 


<P 


$ 


d R $> 


a z $ 


d^ 







66 


66 


64 


1.0 x 10~ 2 


3.4 x 10~ 3 


2.7 x 10~ 3 


4.0 x 10" 


-18 




66 


66 


128 


1.0 x 10~ 2 


3.4 x 10~ 3 


2.7 x 10~ 3 


4.0 x 10" 


-18 




130 


130 


128 


34 x 10~ 3 


1.1 x 10~ 3 


9.7 x 10~ 4 


4.2 x 10" 


-18 




130 


130 


256 


34 x 10~ 3 


1.1 x 10~ 3 


9.7 x 10~ 4 


4.2 x 10" 


-18 


0.5X 


66 


66 


64 


3.3 x 10~ 3 


9.5 x 10~ 4 


7.3 x 10~ 4 


3.3 x 10" 


-4 




66 


66 


128 


2.6 x 10~ 4 


3.6 x 10~ 4 


3.6 x 10" 4 


1.4 x 10" 


-4 




130 


130 


128 


2.3 x 10~ 4 


2.3 x 10~ 4 


1.8 x 10~ 4 


7.0 x 10" 


-5 




130 


130 


256 


9.9 x 10~ 5 


8.6 x 10~ 5 


7.1 x 10~ 5 


3.6 x 10" 


-5 



Table 5 

Maximum Error for Gravitational Potential and Force 

Origin R z <j) $ d R <$> d z <S> <9 $ 

66 66 64 1.3 x 10~ 2 4.7 x 10~ 2 4.8 x 10~ 2 8.3 x 10" 17 

66 66 128 1.3 x 10~ 2 4.7 x 10~ 2 4.8 x 10~ 2 7.7 x 10~ 17 

130 130 128 4.6 x 10~ 3 2.3 x 10~ 2 2.1 x 10~ 2 6.4 x 10~ 17 

130 130 256 4.6 x 10~ 3 2.3 x 10~ 2 2.1 x 10~ 2 7.3 x 10~ 17 

0.5X 66 66 64 8.8 x 10~ 3 1.2 x 10" 1 1.1 x 10" 1 2.4 x 10~ 2 

66 66 128 3.7 x 10~ 3 5.0 x 10~ 2 5.2 x 10" 2 2.2 x 10" 2 

130 130 128 3.4 x 10~ 3 6.0 x 10~ 2 4.7 x 10~ 2 1.3 x 10" 2 

130 130 256 1.0 x 10~ 3 3.3 x 10~ 2 3.2 x 10~ 2 1.8 x 10~ 2 
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For the case where the sphere is centered on the 
origin of the computational grid, the resulting po- 
tential is axisymmetric to machine accuracy. The 
average relative error in the potential and the aver- 
age absolute error in the radial and vertical deriva- 
tives all decrease by a factor of about three as the 
radial and vertical resolutions arc doubled. As ex- 
pected for an axisymmetric mass distribution the 
quality of the solution is independent of the num- 
ber of azimuthal zones. The maximum values of 
the relative error in the potential decrease by a 
factor of about three as well and the maximum 
value in the absolute error of the radial and verti- 
cal derivatives has been cut in half as the number 
of radial and vertical zones doubles. 

When the sphere is placed off axis, the conver- 
gence pattern is much more difficult to recognize. 
For the off-axis test at the highest resolution (the 
same radial and azimuthal resolution that we cur- 
rently use for binary evolutions), we are able to 
obtain a solution that is accurate to one part in 
10 4 , on average, for the potential. Similarly, the 
finite-difference and analytical components of the 
derivatives of the potential agree to better than 4 
decimal places on average. 

5.3. Test of Hydrostatic Equilibrium 

A stringent test of our coupled solution of Pois- 
son's equation and the fluid dynamical equations 
- and one that may seem trivial at first men- 
tion — is how well we are able to maintain hy- 
drostatic equilibrium for a simple system such as 
a spherical polytrope that is placed off axis in the 
grid. While our hydrodynamics implementation 
is conservative with respect to the advection of 
the fluid, there is no guarantee that the total mo- 
mentum is conserved once the action of the La- 
grangian source terms are included. Throughout 
a mass-transfer simulation, the bulk of the fluid 
should remain near hydrostatic equilibrium and 
the correct response of both components to their 
changing mass can be limited by the accuracy to 
which force balance is maintained. 

To perform this test, we have placed a spherical, 
n = 3/2 polytrope of radius R — 0.38 in a cylin- 
drical grid of total radius 1.0, but with a variety 
of different resolutions. The polytrope is centered 
at x « 0.58. In each case, the initial density dis- 
tribution was generated with our SCF code (with 
only one star present and no frame rotation), and 



the initial velocities were zero everywhere. Us- 
ing our full gravitational hydrodynamics code, we 
then permitted the fluid system to evolve in time. 

Over the course of the evolutions, each isolated 
star drifts outwards as if acted on by a constant 
force. This drift is shown in Fig. 6 where we 
have plotted the location of the center of mass 
of the star as a function of time for grids of 
varying resolution. We have normalized the evo- 
lution time to the dynamical time as given by, 
for example, Chandrasckhar (1939). Specifically, 
^dynamical = y/ (3ir) / (16G p) , and for an n = 3/2 
polytrope with central density of unity the average 
density is p — 0.1669. The size and rate of the drift 
decreases as the azimuthal resolution increases. 

As another measure of the quality of the steady 
state equilibrium from these spherical polytropes 
we show a modified virial error, 



VE = 



(W + 311) 

\W\ ■ 



(59) 



in Fig. 7. This differs from the definition given 
in eq. (21) in that we have neglected the kinetic 
energy term, K. As can be seen in Fig. 8 where 
we show the log of \W\, II and K normalized to 
the initial value of \W\ for the highest resolution 
simulation (computed with 130 radial and verti- 
cal zones by 256 azimuthal zones), some peaks 
in kinetic energy, which are noise, are of approxi- 
mately the same size as the sum of W and 311 in 
spite of the fact that the kinetic energy is insignif- 
icant compared to either the thermal or gravita- 
tional energies. Overall, the virial error decreases 
by a factor of approximately 6 from the lowest 
to highest resolution simulation. At the highest 
resolution presented the virial error is 0.05% and 
the polytrope oscillates with amplitude of approx- 
imately 0.02% for 30 dynamical times. This shows 
that the isolated star remains in hydrostatic equi- 
librium to a very high degree of accuracy, even 
when placed off-axis in our computational grid. 

There is no significant improvement in the drift 
of the system center of mass when only the radial 
and vertical resolutions are increased. (In Fig. 6, 
compare the curves for the simulations at reso- 
lutions of 66 radial, 66 vertical by 128 azimuthal 
zones and 130 radial, 130 vertical by 128 azimuthal 
zones.) But there is an improvement in the virial 
error between these two simulations. This sug- 
gests that there are two limiting numerical effects 
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Fig. 5. — Clockwise from the upper left panel, the 
vertical velocity, pressure, the ratio of the entropy 
tracer to the mass density and the mass density are 
shown as a function of position z at time t — 0.247 
for Sod's shock tube problem. The initial discon- 
tinuity was placed at z = —0.1. The computed 
solution for a randomly chosen column of cells at 
constant radius and azimuth is plotted as crosses; 
the solid curves are the analytic solution. 
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Fig. 7. — The virial error, as given by eq. (59) ( the 
virial error with the kinetic energy term omitted) 
is plotted as a function of the number of dynamical 
times for the off-axis spherical polytrope described 
in §5.3. The meaning of the different curves is the 
same as in Fig. 6. 
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Fig. 6. — Distance from origin to the center of 
mass as a function of time (measured in dynami- 
cal times) for the off-axis spherical polytrope de- 
scribed in §5.3. The different curves correspond 
to calculations performed at the indicated resolu- 
tion, in terms of the number of radial by vertical 
by azimuthal zones. The curves representing the 
simulations at resolutions of 66 x 66 x 128 and 
130 x 130 x 128 lie on top of one another. 
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Fig. 8. — The logarithm of the three components 
to the virial error as given in eqs. (22)— (24) nor- 
malized to the initial value of the gravitational po- 
tential energy is plotted as a function of time for 
the off-axis spherical polytrope evolved in a grid 
of resolution 130 x 130 x 256. 
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in play. One dictates the resolution of the equilib- 
rium state itself; the other causes a displacement 
of that equilibrium state. The former effect con- 
verges with the finite-difference size isotropically, 
while the latter depends only on the azimuthal 
resolution. When trying to resolve a highly non- 
axisymmetric object, such as an off-axis sphere, 
within a uniform cylindrical coordinate grid, dif- 
ferent parts of the star are resolved to varying de- 
grees and it is not surprising that the convergence 
of the numerical solution is not describable in sim- 
ple terms. 

6. Benchmark Simulations 

In this section we present results from two sim- 
ulations of detached binaries that we have per- 
formed to ascertain the precision with which we 
can expect to carry out future simulations of 
semi-detached binary systems (systems undergo- 
ing mass transfer). One binary is an equal mass 
system with identical components (see Model 5 in 
Table 1 and Figs. 2-3; hereafter referred to as the 
EB system) and the other system has a mass ratio 
q = 0.8436 (see Model 6 in Table 1 and Figs. 2- 3; 
hereafter referred to as unequal binary or UB sys- 
tem) . The EB system was constructed to resemble 
the single star used for the test of hydrostatic equi- 
librium in §5.3. This enables us to compare the 
systematic errors in the case of a binary system 
given the errors observed when only gravity, pres- 
sure and the curvature force came into play. Each 
component of the EB system differs from the iso- 
lated, spherical star in that each is flattened by the 
synchronous rotation of the system and tidally dis- 
torted by its companion, but the components have 
a comparable size, in terms of grid cells, and the 
same central density and polytropic index as the 
isolated sphere. 

Previous simulations of equal-mass barotropic 
stars have shown that it is important to conduct 
the evolutions in a frame of reference that renders 
the binary as close to static as possible in order 
to minimize the effects of numerical diffusion aris- 
ing from the finite accuracy of Eulerian advection 
schemes (New & Tohline 1997; SWC). With this 
in mind, our EB and UB simulations have been 
conducted in a frame of reference rotating with 
the orbital angular velocity £1 of the system, as 
obtained by our SCF technique. 



In dealing with unequal-mass systems we have 
discovered another subtle, but important issue 
that should be addressed with care when "trans- 
porting" an initial hydrostatic model from the grid 
of the SCF code into the grid of the hydrodynam- 
ics code. During each SCF iteration, the system's 
center of mass is not fixed to any location beyond 
the fact that, by symmetry, it must lie along the 
line of centers. In general, then, we must trans- 
late the density field as we introduce it into the 
hydrocode so that the system center of mass coin- 
cides with the z-axis, which is taken to be the ro- 
tation axis for the hydrodynamic evolution. If we 
could perform this translation perfectly, all initial 
fluid velocities would be identically zero relative to 
the hydrodynamic reference frame. Because of the 
inherent symmetry of an equal-mass binary sys- 
tem, this was in fact the case for our EB system 
by construction. For the UB system, however, the 
center of mass of our converged SCF model was 
displaced by a small distance from the rotation 
axis. Specifically, R CO m = 2.842 x 10~ 6 x, which 
corresponded to only 4x 10 _4 Ai?, where AR is the 
radial extent of each grid cell. As we introduced 
the SCF model into the hydrodynamical grid, we 
therefore also ascribed nonzero velocities as initial 
conditions according to the relation, 



n x R r 



(60) 



Because the displacement R CO m was quite small 
for our UB system, the initial velocities prescribed 
through eq. (60) were also very small. Neverthe- 
less, it was necessary to include them in order to 
achieve the best possible steady-state configura- 
tions corresponding to the stars following circular 
orbits. This implies a uniform initial velocity for 
the system (see further discussion below). 

After the binary models were introduced into 
the hydrodynamics code, both were evolved 
through more than five orbits. (See the first row 
of Table 6, where the total evolution time for both 
simulations is tabulated in units of each system's 
orbital period P.) As is recorded in the last three 
rows of Table 6, the EB system was run on 64 
nodes of a Cray T3E 600 for a total of 173 wall- 
clock hours (that is, the simulation required on 
average 2409 processor-hours per orbit) and the 
UB system was run on 8 dual processor nodes of 
an IBM SP 3 for a total of 265 wall-clock hours 
(that is, the simulation required on average 819 
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processor- hours per orbit). Many different diag- 
nostic parameters were followed throughout both 
evolutions in order to assess the quality of the ini- 
tial SCF models and to determine with what ac- 
curacy the hydrodynamical equations were being 
integrated forward in time. In the following para- 
graphs, we present the time-evolutionary behavior 
of a number of these key physical parameters. 

6.1. Stars in Hydrostatic Balance 

Throughout both evolutions, the individual 
stellar components were largely static and re- 
mained well within their respective Roche lobes. 
In an effort to illustrate this, Fig. 9 shows as a 
function of time the computed Roche lobe volume 
(dashed curve) and the volumes (solid curves) oc- 
cupied by material more dense than 10 _1 , 10 -2 , 
10~ 3 , 10~ 4 , and 1CP 5 for one component of the 
EB system. (For reference, the initial SCF density 
fields have values of a few times 10~ 5 at the edge 
of the stars; see the isodensity contours drawn in 
Figs. 2 and 3.) The same information is plotted 
in Figs. 10 and 11 for the secondary and pri- 
mary components, respectively, of the UB sys- 
tem. These figures illustrate that the rotationally 
flattened and tidally distorted models generated 
by our SCF code exhibit excellent detailed force 
balance throughout their three-dimensional struc- 
tures, and that there is an excellent match be- 
tween the algorithmic expressions that determine 
an equilibrium state in the SCF code and force 
balance in the hydrodynamics code. 

6.2. Mass Conservation 

In an effort to determine how well mass is con- 
served throughout an evolution for each star, indi- 
vidually, as well as for the system as a whole, we 
tracked three separate volume integrals over the 
mass density: M\, defined as the mass bound to 
the primary; M2, defined as the mass bound to the 
secondary; and M onvo i opo , defined as the mass that 
lies outside of both stars but inside the bound- 
aries of the computational grid. As is illustrated 
by frames 5 and 6 of Figs. 2 and 3, in the initial 
state it is easy to evaluate these three integrals be- 
cause the edges of the two stars are well-defined. 
Specifically, when normalized to each system's to- 
tal mass, Mi = (1 + g) _1 , M 2 = q(l + g) _1 , and 
A/e nve iope = 0, where q is the system mass ra- 
tio given it Table 1. But because the stars are 
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t/P 



Fig. 9. — The Roche volume (dashed curve) and 
volume occupied by material more dense than 
KT 1 , 1CT 2 , 1CT 3 , 1(T 4 and 1(T 5 (solid curves 
from bottom to top) are plotted as a function of 
the orbital time for one component of the EB sys- 
tem. 




Fig. 10. — The Roche volume (dashed curve) 
and volume occupied by material more dense than 
10 -1 , 10~ 2 , 10~ 3 , 10~ 4 and 10~ 5 (solid curves 
from bottom to top) as a function of the orbital 
time for the secondary component of the UB sys- 
tem. 
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Table 6 
Quantities of Interest for Benchmark Simulations 

Quantity Equal Mass Binary (EB) Unequal Mass Binary (UB) 



AMi 

M 
AM 2 

ik 
. M 

( — ) 

V a / secular 

— 

V a / epic 

Jz 

Machine 
Processors 

^WallClock 



' cpicyclic 



5.314 

-9.0 x 10~ 6 

-1.0 x 10~ 5 

-1.9 x 10- 5 

-2.9 x 10~ 4 

5.0 x 10~ 4 

+1.1 x 10~ 4 

Cray T3E 600 

64 

173 hours 



5.178 
-3.0 x 10- 5 
-1.1 x 10~ 6 
-1.4 x 10~ 5 
-1.9 x 10~ 4 

2.2 x 10~ 4 
+ 1.5 x 10~ 4 

IBM SP3 
16 

265 hours 



S 0.15 

01 

E 
-± 

> 0.10 
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Fig. 11. — The Roche volume (dashed curve) 
and volume occupied by material more dense than 
10~ 5 , 10~ 4 , 10~ 3 , 10~ 2 and lO" 1 (solid curves 
from bottom to top) as a function of the orbital 
time for the primary component of the UB system. 



being modeled on a discrete computational mesh 
that does not conform precisely to their shape, and 
because the acceleration of each fluid element in 
the computational mesh is being determined by 
finite-difference (rather than continuous differen- 
tial) representations of gradients in the pressure 
and gravitational fields, as each system evolves hy- 
drodynamically the surfaces of the stars become 
less sharply defined and some spreading of mate- 
rial inevitably occurs. (In these benchmark evo- 
lutions, this is evidenced for example by the very 
small, but finite oscillations in the "volumes" oc- 
cupied by the stars that are displayed in Figs. 9 
- 11.) In practice, then, during each evolution we 
determine whether material in each grid cell be- 
longs to cither star or the "envelope" by compar- 
ing the binding energy of the fluid in each cell to 
the average binding energy of the layer of cells at 
the surface of each star. In this context, we de- 
fine the surface of each star to be the layer of cells 
where the mass density falls below 10~ 5 in our 
normalized units, which corresponds to the lowest 
density level attained in the initial SCF models. 
The mass of the envelope is dominated by material 
from the surface of the stars even though there is 
a minimum "vacuum" density level of 1.0 x 10~ 10 
enforced by the code to maintain numerical sta- 
bility. The total mass of the vacuum material is 
over a million times smaller than the mass of ei- 
ther stellar component and does not impact the 
physics of these simulations. 
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Four curves are drawn (each at two quite differ- 
ent scales) in Figs. 12 and 13 to document how well 
mass is conserved in the EB and UB simulations, 
respectively. The masses have all been normalized 
to the total system mass, so in Fig. 12 the mass 
of both the primary (Mi) and the secondary (M2) 
stars is initially exactly 0.5; the total binary mass 
is initially exactly 1; and the "envelope" mass is 
initially exactly 0. In Fig. 13, the total mass and 
the "envelope" mass also are initially 1 and 0, re- 
spectively, but the mass of the primary initially 
is (1 + q)^ 1 = 0.5424 and the mass of the sec- 
ondary initially is q(l + q)^ 1 = 0.4576. Plotted 
on a normal, linear scale, all four of these curves 
are perfectly flat in both figures. This demon- 
strates that the mass of both stars, as well as the 
aggregate system mass, is conserved to very high 
precision throughout the EB and UB simulations. 
Again, this is evidence that the initial models were 
in excellent detailed force balance and the hydro- 
dynamics code is evolving the systems forward in 
time in a physically realistic manner. 

Although mass is conserved to very high preci- 
sion, it is not absolutely constant throughout the 
evolutions. In the four insets of Figs. 12 and 13, we 
have magnified the vertical mass scale by roughly 
four orders of magnitude in order to show that 
there is a very tiny, but measurable, secular de- 
crease in the total system mass and in the mass of 
both stellar components over the course of the sim- 
ulations. In each inset, we plot the relevant mass 
minus its value in the initial state (time t = 0), 
normalized to the total system mass. These inset 
plots show that the system mass decreases by ap- 
proximately one part in 10 4 over five orbits — that 
is, about 0.002% per orbit — with the mass loss 
from each star accounting for roughly half this to- 
tal. In rows 2 — 4 of Table 6 we have recorded for 
both evolutions more precise values of the frac- 
tional mass that is lost, on average, each orbit 
from the primary (Mi), the secondary (M2), and 
the system as a whole (M). We have determined 
that this mass is very slowly lost as a result of 
the development of a small, but nonzero flow of 
low-density material off of the stars, through the 
envelope and, ultimately, off of the grid. (After an 
initial drop, the envelope mass remains approxi- 
mately constant, suggesting that this outward flow 
has settled into a nearly steady state.) As is dis- 
cussed more fully in §7, below, this small but de- 
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Fig. 12. — Masses, normalized to the total system 
mass, plotted as a function of time, in units of the 
orbital period, for the EB simulation. Top curve: 
Total binary mass (M). Middle two curves: Mass 
of the primary (Mi) and secondary (M2) stars. 
Bottom curve (essentially at zero): Mass of the 
"envelope," as defined in §6.2. Inset plots show the 
difference between the indicated mass component 
and its initial value in units of the initial total 
mass. 
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tectable rate of mass loss from detached equilib- 
rium binaries imposes a straightforward limit on 
the mass-transfer rates that we will be able to re- 
liably model in future simulations that involve dy- 
namical mass transfer. 
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Fig. 13. — Masses, normalized to the total system 
mass, plotted as a function of time, in units of the 
orbital period, for the UB simulation. Top curve: 
Total binary mass (M). Middle two curves: Mass 
of the primary (Mi) and secondary (M2) stars. 
Bottom curve (essentially at zero): Mass of the 
"envelope," as defined in §6.2. Inset plots show the 
difference between the indicated mass component 
and its initial value in units of the initial total 
mass. 



6.3. Minimal Center-of-Mass Motion 

During both simulations we also tracked as a 
function of time the position of the center of mass 
of each binary component and the position of the 
center of mass of the system as a whole. The 
equatorial-plane trajectories of these three cen- 
ters of mass for the EB and UB evolutions are 
shown, respectively, in Figs. 14 and 15, as viewed 
from our computational reference frame — that 
is, from a frame rotating with the orbital angular 
velocity of the system, as determined for the ini- 
tial state by the SCF code. In the uppermost plot 
of each figure, which has been drawn at a scale 
(—1 < x < +1; — 1 < y < +1) to include the en- 
tire mass of the system, the three separate center 
of mass trajectories appear to be small dots with 
little or no discernible structure. (When plotted 
in the inertial reference frame on this scale the 
trajectories of the two stars are indistinguishable 
from circles.) This illustrates that, even after five 
orbits, the centers of mass of the two stars and of 
the system as a whole essentially have not moved 
from their initial positions. This provides addi- 
tional strong confirmation that our SCF code pro- 
duces excellent initial states and that the hydro- 
dynamical equations are being integrated forward 
in time in a physically realistic manner. 

In the bottom three plots of Figs. 14 and 15, we 
have magnified a small region around each of the 
center of mass trajectories — expanding the lin- 
ear scale of the uppermost plot in each figure by a 
approximately a factor of 15 and 45, respectively 
These magnified views reveal that, although it is 
very small, there is measurable motion of the cen- 
ters of mass in both evolutions. In the bottom, 
lcfthand plot we also have shown the size of our 
radial grid spacing, AR = 7.87 x 10~ 3 . This in- 
dicates the characteristic size of our discretization 
and emphasizes how small the motion of each cen- 
ter of mass is. In the UB evolution, for example, 
the motion of all three centers has been confined 
within a single grid cell through five full orbits. 
Furthermore, the smooth spiral trajectory of the 
UB system center of mass (bottom, middle plot in 
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Fig. 14. — From the EB simulation, equatorial- 
plane trajectories are plotted for the center of mass 
of the system and the centers of mass of both stel- 
lar components through just over 5 orbits in the 
corotating frame of reference. Insets (from left to 
right) show magnified views of the trajectories for 
the secondary star, the system as a whole, and 
the primary star. We have subtracted off the ini- 
tial coordinates for the inset plots and have, for 
reference, indicated the size of one grid cell. 



Fig. 15. — From the UB simulation, equatorial- 
plane trajectories are plotted for the center of mass 
of the system and the centers of mass of both stel- 
lar components through just over 5 orbits in the 
corotating frame of reference. Insets (from left to 
right) show magnified views of the trajectories for 
the secondary star, the system as a whole, and 
the primary star. We have subtracted off the ini- 
tial coordinates for the inset plots and have, for 
reference, indicated the size of one grid cell. 



2.1 



Fig. 15) has an understandable, physical origin. 
As viewed in the inertial frame, this particular tra- 
jectory appears as a straight line whose direction 
and magnitude is consistent with the overall sys- 
tem velocity prescribed as initial conditions from 
eq. (60). In the EB evolution, due to the symmetry 
of the initial model, the drift of the system center 
of mass is extremely small, remaining unnotica- 
ble even on the magnified plot. In the magnified 
plots, the trajectory of the center of mass of each 
individual star shows both a gradual drift in the 
y-direction, and a small oscillatory motion in the 
x-direction. The vertical drift is mostly an indica- 
tion that the binary's actual orbital frequency is 
slightly different from the value (given by the SCF 
code) that we used for the rotation frequency of 
the computational grid. The oscillations in x rep- 
resent epicyclic motion and indicate that the bi- 
nary orbit is not precisely circular. Since both the 
drift and the epicyclic motion can be understood 
in physical terms, their small amplitudes tell us 
more about the quality of the initial model than 
about the limiting accuracy of our finite-difference 
scheme. 

Unlike in the single star case presented in §5.3, 
there is no evidence of a systematic outwards force 
on either star in the UB or EB systems, despite the 
fact that the systems have evolved for the equiv- 
alent of approximately 90 dynamical times. It 
appears as though the introduction of a rotating 
frame of reference and the associated centrifugal 
potential and Coriolis force has provided a feed- 
back mechanism that acts to limit the systematic 
imbalance discussed previously. 
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6.4. Binary Separation 

A plot of the binary separation a as a func- 
tion of time, as shown in Figs. 16 and 17 for the 
EB and UB evolutions, respectively, provides an- 
other way to assess the global behavior of these 
systems. Here, the separation a is defined as the 
distance between the centers of mass of the two 
stars. Notice that, on a linear scale that extends 
from to 1 (in units normalized to the each sys- 
tem's initial separation), the plot of a(t) is indis- 
tinguishable from a perfectly horizontal line. This 
illustrates that, to a very high degree of accuracy, 
these benchmark simulations of detached binary 
systems produce stable, circular orbits. 

Again, though, if we examine these plots in 



Fig. 16. — The orbital separation, normalized to 
its initial value, as a function of orbital time for 
the EB system. 
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Fig. 17. — The orbital separation, normalized to 
its initial value, as a function of orbital time for 
the UB system. 



finer detail, we see that both evolutions exhibit 
a very small but quantifiable departure from per- 
fect circular orbital motion. For example, in the 
insets to Figs. 16 and 17, we have reploted a(t) 
with the vertical scale magnified by roughly a fac- 
tor of 400. These insets show that in both evo- 
lutions there is a very slow, secular decrease in 
the orbital separation and, in addition, a(t) dis- 
plays low-amplitude oscillations having a period 
approximately equal to one orbital period. The 
oscillations in a arise from the same epicyclic mo- 
tion that was seen in the plots (Figs. 14 and 15) of 
the center of mass motion of the individual stars, 
but the amplitude of this motion is easier to mea- 
sure here. In units of the initial orbital separa- 
tion, the EB system has an epicyclic amplitude 
(Aa/a)epi cyc ii c w 5 x 10~ 4 ; the UB system exhibits 
an epicyclic amplitude about half this size. The 
slow, secular decay of the orbits occurs at a rate 
(Aa/a) secu i ar ~ 2.9 x 10~ 4 per orbit in the EB 



system, and at a rate Aa/a w 1.9 x 10 per or- 
bit in the UB system. These orbital decay rates 
and epicyclic amplitudes have been recorded in the 
fifth and sixth rows of Table 6. 

6.5. Angular Momentum Conservation 

Finally, in Figs. 18 and 19 we show the behav- 
ior as a function of time of the z-component of 
each system's total angular momentum. As was 
true with our plots of the orbital separation, on 
a linear scale that extends from to 1 (in units 
normalized to the each system's initial total an- 
gular momentum), the plot of J z (t) is indistin- 
guishable from a perfectly horizontal line. This 
illustrates that these benchmark simulations glob- 
ally conserve angular momentum to a very high 
degree of accuracy. When we magnify the verti- 
cal scale by approximately a factor of 1000, as has 
been done to produce the insets to Figs. 18 and 
19, we see that angular momentum is not, in fact, 
perfectly conserved. Evidently both systems gain 
angular momentum at a very slow rate; in the EB 
simulation, AJ Z /J Z « 1.1 x 10 -4 per orbit, and in 
the UB simulation AJ Z /J Z « 1.5 x 10~ 4 per or- 
bit. These rates have been recorded in the seventh 
row of Table 6 and will be referred to again in §7, 
below, when we summarize the limiting accuracy 
with which we expect to be able to model physical 
mass-transfer events using our simulation tools. 
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Fig. 18. — The z component of total angular mo- 
mentum, normalized to its initial value, as a func- 
tion of time for the EB system. 



Fig. 19. — The z component of total angular mo- 
mentum, normalized to its initial value, as a func- 
tion of time for the UB system. 
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6.6. Overview 



We should emphasize that the hydrodynamics 
code as described in §4 and utilized in these bench- 
mark simulations has evolved through many stages 
from the version of the code that was used sev- 
eral years ago by New & Tohline (1997) to inves- 
tigate the equal-mass, binary merger problem. A 
number of improvements were made in the code 
in order to bring it to its present level of perfor- 
mance. Figure 20 is presented here in an effort 
to illustrate how certain key modifications in the 
code affected its performance. Each row of frames 
in this figure shows results from an evolution of 
the same unequal-mass (UB) binary system that 
was used in our benchmark simulation, but as pro- 
duced by four separate versions of the code. The 
curves drawn in the four frames on the left-hand 
side of Fig. 20 show the same type of information 
as has been presented in Fig. 10 for the bench- 
mark UB evolution: Four separate volumes for 
the secondary star (solid curves) and its Roche 
volume (dashed curve) are plotted as a function 
of time, in units of the orbital period. The four 
frames on the right-hand side of Fig. 20 show the 
same type of information as has been presented 
in the top-most frame of Fig. 15: Center-of-mass 
trajectories of the two stars and of the system as 
a whole, as viewed from the rotating frame of ref- 
erence. (The bottom-most frames are taken from 
the benchmark UB simulation and therefore are 
drawn directly from Figs. 10 and 15.) 

The results shown in the top-most frames of 
Fig. 20 come from an early version of the code 
in which we replaced the gradient of the pressure 
with the gradient of the enthalpy. This ensured 
that the initial structure of each star, as deter- 
mined by the SCF code, was in good force bal- 
ance after being introduced into the hydrocode. 
However, as the two figures from this evolution il- 
lustrate, we still observed a slow expansion of the 
secondary star; the orbit itself developed a signif- 
icant epicyclic amplitude; and after about three 
orbits, the Roche lobe was encroaching on the sur- 
face of the secondary. The second row of frames 
comes from a simulation in which the number of 
azimuthal zones was doubled — from 128 to 256 
zones over the full 27r radians. This modification 
improved somewhat the mean motion of the cen- 
ters of mass (although it did not significantly re- 
duce the amplitude of the epicyclic motion) . Most 
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Fig. 20. — Each of the four pairs of plots shown 
here has been derived from the UB simulation as 
modeled with one of four separate versions of our 
hydrodynamic code. In each pair, the plot on 
the left is directly analogous to Fig. 10, showing 
as a function of time the Roche volume (dashed 
curve) and four interior volumes (solid curves) for 
the secondary star; the plot on the right is di- 
rectly analogous to the unmagnificd plot in Fig. 
15, showing the trajectories of the center of mass 
of the secondary (left), system as a whole (center), 
and primary (right). The pairs of plots are shown 
chronologically from the top to the bottom, with 
significant improvements in the code being made 
between each recorded UB simulation. See §6.6 for 
a description of these various code improvements. 
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significantly, however, doubling the angular grid 
size improved the resolution and, hence, the de- 
termination of force balance in each star. As a 
result, expansion of the secondary star was no- 
ticeably reduced. The third row of frames shows 
that motion of the centers of mass was drasti- 
cally reduced when we modified our algorithm to 
make the integration scheme more properly time- 
centered. This change did not noticeably reduce 
the rate of expansion of the secondary, but it did 
significantly reduce the amplitude of oscillations 
in the Roche lobe volume. Finally, by introducing 
artificial viscosity into the equations of motion in 
order to mediate the weak shocks at the surface 
of the stars (which also involved a re-centering of 
all momentum densities to the cell locations spec- 
ified in Table 3), the entire structure of both stars 
became much more robust. In particular, as the 
left-hand frame of the last row shows (see also Fig. 
10), this code modification completely eliminated 
the short timescale wiggles in the volumes of the 
secondary; overall expansion of the secondary also 
was reduced to an imperceptible level. Simulta- 
neously, for the first time, we ascribed a small 
nonzero velocity to the initial state as given by 
cq. (60). This change further reduced the motion 
of the centers of mass — to the level illustrated by 
Fig. 15. 

It is reasonable to ask whether the three prin- 
cipal spurious effects that remain in our bench- 
mark simulations — the slow decay of the orbits, 
the slow gain of angular momentum, and the slow 
loss of mass from the stars — are at least mutu- 
ally consistent on physical grounds. For centrally 
condensed binaries, a point mass approximation 
(the Roche approximation) is usually sufficient to 
discuss the orbital evolution. This approximation 
predicts a simple relationship between the changes 
of mass, angular momentum and binary separa- 
tion, namely 
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where J com = M\M2\/GajM is the total center- 
of-mass angular momentum in the point mass ap- 
proximation. Numerical values for each of the 
terms in this expression can be obtained from the 
data shown in Table 6. In particular, we sec that in 
both benchmark simulations the three mass terms 
approximately cancel each other out. But while 



the magnitude of (Aa/a) socu i ar is roughly the same 
as the magnitude of AJ Z /J Z , their signs are differ- 
ent. That is, the angular momentum of the system 
is slowly increasing while the binary separation 
is slowly decreasing. This is clearly inconsistent 
with the expectations of the simple Roche model. 
A more accurate expression for the total angular 
momentum of the binary would be 



f-'bin — ^ci 



JifiiH- J 2 ft 2 



(62) 



where Ii and Oj are the moments of inertia and in- 
ertial frame angular velocities of the binary com- 
ponents, assuming they all rotate around the same 
z-axis. Even if one takes into account the con- 
tributions of spin angular momenta, the changes 
observed remain inconsistent and must therefore 
be attributed to spurious numerical effects at a 
level of 10 -4 per orbit arising from the inevitable 
error terms present in our finite-difference repre- 
sentation of the fluid equations. What we have 
attempted to do here is quantitatively document 
the magnitude of these numerical effects in the 
highest practical resolution possible at the present 
day for simulations of detached binaries where the 
character of the ideal solution is well understood 
beforehand. Furthermore, we can not accurately 
predict the evolution of mass transferring binaries 
where the mass transfer rate, AM /M per orbit is 
< few x 10~ 5 using simulations at the resolution 
presented here. There are however a wide variety 
of systems (the initial mass transfer event in an Al- 
gol progenitor, or the onset of common envelope 
evolution in the progenitors of many types of bi- 
naries, or the formation of Type la supcrnovae for 
example) that are expected to exceed our thresh- 
old resolution limit for mass transfer. At a suffi- 
ciently high mass-transfer rate, the mass transfer 
itself will drive the evolution of the orbital param- 
eters and Roche geometry at a rate higher than 
the numerical limits demonstrated here. 

7. Conclusions 

In this paper we have presented a practical SCF 
algorithm for constructing self-consistent poly- 
tropic binaries with unequal masses that satisfy 
the condition of hydrostatic equilibrium to a high 
degree of accuracy. This three-dimensional SCF 
algorithm is based largely on the technique first 
described by Hachisu (1986), but to our knowledge 
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this is the first time the technique has been used 
to generate unequal mass binary systems as input 
for a hydrodynamical simulation. Our two bench- 
mark simulations (described in §6) clearly indicate 
that this SCF algorithm can provide superb ini- 
tial states for investigations into the dynamical 
stability of close binary systems. We emphasize 
that, in addition to generating models of close bi- 
nary systems that are detached — like the EB and 
UB systems constructed for our two benchmark 
simulations — as shown above in Figs. 2 and 3 
this technique also can be used to generate close 
binary systems that are semi-detached or, in the 
limit of identical components, in contact. 

We have also detailed our gravitational hydro- 
dynamics code and presented results from key 
tests of the stability and accuracy of the hydrody- 
namics algorithm, the solution of Poisson's equa- 
tion, and the coupled solution required to evolve 
an isolated, spherical polytropc that is placed off- 
axis in a cylindrical grid. From these test cases 
it is apparent that a number of subtle numer- 
ical issues arise when a highly nonaxisymmct- 
ric body is evolved via an explicit integration of 
finite-difference equations on a cylindrical compu- 
tational grid. It also appears however, that these 
effects can be made manageably small by increas- 
ing the resolution used to treat the system of in- 
terest. 

We have evolved two detached binary systems 
(one with a mass ratio q = 1, the other with a 
mass ratio q — 0.8436) through more than five 
orbits in order to benchmark in detail the capabil- 
ities of our simulation tools. Even though the indi- 
vidual stellar components generated by our SCF 
code are significantly rotationally flattened (due 
to the synchronous rotation of the initial states) 
and tidally distorted (by their close binary com- 
panion), these benchmark simulations show that 
the stars are in almost perfect hydrostatic equi- 
librium. Throughout each binary evolution, our 
hydrodynamics code conserves mass and angular 
momentum to a very high degree of precision; as 
viewed from a frame of reference that is rotating 
with the initial orbital frequency of the binary, the 
centers of mass of the two stellar components and 
of the system as a whole remain virtually station- 
ary; and a plot of the binary separation as a func- 
tion of time shows that the stellar orbits are almost 
indistinguishable from circles. This gives us con- 



siderable confidence that these numerical tools can 
be used to examine the stability of close binary 
systems against both tidal and mass-transfer in- 
stabilities, and to begin to accurately model mass 
transfer in semi-detached systems. 

As has been summarized in Table 6, from our 
benchmark evolutions we have been able to deter- 
mine in quantitative terms the level of accuracy 
with which our hydrodynamical code conserves 
mass, conserves angular momentum, and is able 
to represent and maintain a circular binary orbit. 
Mass is conserved to roughly 0.002% per orbit; 
angular momentum is conserved to a level of 0.01 
- 0.02% per orbit; the binary separation remains 
constant to a few parts in 10 4 per orbit; and each 
orbit exhibits an epicyclic amplitude (measured 
relative to the orbital separation) of .02-. 05%. 

We are unaware of any other group that is at- 
tempting to study the onset of mass-transfer in- 
stabilities in unequal-mass binaries with a gravi- 
tational hydrodynamics code, like ours, that fully 
resolves both stellar components. Hence, there are 
no published numbers against which to compare 
ours for the UB evolution. However, we can fairly 
compare the results of our EB evolution against 
the recent study published by SWC of equal-mass 
close binary systems. Their Fig. 10 illustrates 
that, after following one stable binary system 
through approximately 6 orbits (we assume, based 
on the SWC discussion, that P = 1.7 — 2 ms), they 
have been able to conserve angular momentum to 
a level of about 0.2% per orbit. And their Fig. 
14 shows four stable orbits with epicyclic ampli- 
tudes (measured relative to the binary separation) 
~ 0.3 — 1%. We conclude that, at least in these 
two respects, our simulations improve on the SWC 
models by roughly one order of magnitude. SWC 
do not comment on their level of mass conserva- 
tion; and, because of the visible epicyclic motions 
in their Fig. 14, it is difficult to ascertain to what 
degree the binary separation either decreases or 
increases with time over the course of their evo- 
lutions. In both of our evolutions, however, the 
center-of-mass motion of our stars appears to be 
significantly less than the center-of-mass motion 
depicted for SWC's preferred integration scheme 
in the top-left panel of their Fig. 7. 

The small, but measurable changes in mass, 
angular momentum, and binary separation doc- 
umented here in Table 6 set limits on the types of 
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mass-transfer events that we will be able to model 
with confidence using our present simulation tools. 
For example, if we were to try to model an instabil- 
ity that leads to a flow through the binary's LI La- 
grange point with a mass-transfer rate lower than 
one part in 10 6 per orbit, the physical exchange of 
material between the binary components would be 
swamped by the unphysical process that is caus- 
ing our stars to lose mass to the "envelope" at a 
rate of one to two parts in 10 5 per orbit. If the 
depth of contact between the Roche lobe and the 
surface of the donor star is not sufficient, epicyclic 
motion in the orbit will tend to shut off the mass- 
transfer during part of each orbit. Also, a drift 
in the center of mass of the system can impose a 
limit on the length of time that the binary can be 
evolved before the motion of the binary through 
the grid becomes problematic. We will have to 
contend with all of these issues as we move to 
the next level of our investigation and introduce 
a semi-detached system from our SCF code into 
our hydrodynamical code. We expect nevertheless 
to find a wide range of interesting binary systems 
whose dynamical evolution can be simulated in a 
fully self-consistent fashion through a reasonably 
large number of orbits using the tools that have 
been described in this paper. 

As we have documented in Table 6, the calcu- 
lation of one orbit takes about 33 wall-clock hours 
when utilizing 64 processors of the Cray T3E 600, 
and using 16 processors of the newer IBM SP-3, 
the calculation of one orbit takes about 51 hours. 
The computational workload of a mass-transfer 
simulation is therefore within the reach of cur- 
rent, state of the art, parallel computers given the 
linear scaling of our gravitational hydrodynamics 
code with the number of processors even at a res- 
olution greater than presented here. We note that 
the amount of work performed can be reduced sig- 
nificantly if need be by, for example, freezing the 
gravitational potential until the mass distribution 
has changed significantly as done by Stone & Nor- 
man (1992). The solution of Poisson's equation 
represents about a quarter of the total execution 
time. 

We have been able to estimate the mass trans- 
fer rate required to bring the simulation above the 
level of the noise observed in our benchmark sim- 
ulations. We find that this value is ~ few x 10~ 5 
of the donor's initial mass over an orbital period. 



As discussed in §2, the mass transfer rate should 
scale as a high power of the degree of over-contact 
(as the cube for an n = 3/2 polytropc); further- 
more, for a case where q > 1, that is, the donor 
is initially the more massive star, the degree of 
over-contact will naturally increase with time as 
the star expands and its Roche lobe shrinks. Pro- 
vided that such a binary system can begin mass- 
transfer, the amplitude of the mass transfer rate 
should inevitably reach values higher than indi- 
cated above. Since motion of the center of mass 
of the binary system has been confined to a region 
well within a single computational grid cell even 
after 5 orbits, we are confident that future evo- 
lutions can be followed with confidence through 
at least 20 — 30 orbits, given sufficient comput- 
ing resources. As discussed in the introduction of 
this paper, we understand that the mass transfer 
rates considered here are much larger than those 
found in what are considered typical examples of 
interacting binaries. The methods presented here 
are not applicable to the stable mass transfer ob- 
served in cataclysmic variables or other long lived 
systems, but should serve very well to investigate 
stages of evolution of their progenitors and tran- 
sient events such as the onset of dynamical mass 
transfer and its stability. 
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